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ABSTRACT 

We investigate the long-term evolution of the inclinations of the known classi- 
cal and resonant Kuiper belt objects (KBOs). This is partially motivated by the 
observed bimodal inclination distribution and by the putative physical differences 
between the low- and high-inclination populations. We find that some classical 
KBOs undergo large changes in inclination over gigayear timescales, which means 
that a current member of the low-inclination population may have been in the 
high-inclination population in the past, and vice versa. The dynamical mecha- 
nisms responsible for the time- variability of inclinations are predominantly dis- 
tant encounters with Neptune and chaotic diffusion near the boundaries of mean 
motion resonances. We reassess the correlations between inclination and physical 
properties including inclination time- variability. We find that the size-inclination 
and color-inclination correlations are less statistically significant than previously 
reported (mostly due to the increased size of the data set since previous works 
with some contribution from inclination variability). The time- variability of in- 
clinations does not change the previous finding that binary classical KBOs have 
lower inclinations than non-binary objects. Our study of resonant objects in the 
classical Kuiper belt region includes objects in the 3:2, 7:4, 2:1, and eight higher- 
order mean motion resonances. We find that these objects (some of which were 
previously classified as non- resonant) undergo larger changes in inclination com- 
pared to the non-resonant population, indicating that their current inclinations 
are not generally representative of their original inclinations. They are also less 
stable on gigayear timescales. 

Subject headings: Kuiper belt: general, planets and satellites: djTiamical evolu- 
tion and stability, chaos 



1. Introduction 



The orbital structure of the Kuiper belt contains information about both the formation 
of the solar system and the Kuiper belt's subsequent evolution under the current solar system 
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architecture. The effect of the latter on the Kuiper belt must be well understood in order 
to determine how planetary formation and mi gration sculpted the orb ital distribution of the 
primordial Kuiper belt (see review article by iMorbidelli et al.l (j2008[ )). Understanding the 
evolution of the inclination distribution of the classical Kuiper belt i s of particular interest 
because the bimodal inclination distribution found by iBrownl (120011 ) is an indication that 
classical Kuiper belt objects (CKBOs) might actually be comprised of two overlapping popu- 
lations: a primordial, low-incl ination population, and a dynamically excited, high-inclination 
population (see, for example, iKavelaars et al.l ( l2008l )). 



There have been se veral studies attempt ing to identify an explanation for the bimodal 
inclination distribution. iKuchner et al.l ( 120021 ) investigated the inclination-dependence of the 
long-term dynamical stability of Kuiper belt objects by performing numerical integrations 
of test particles with semimajor axes spread from 41 to 47 AU with uniform inclination, 
up to ~ 30°. They found that in the inner part of the Kuiper belt, secular resonances 
systematically destabilize low-i test particles, resulting in an inclination distribution skewed 
toward higher z; however, they did not find any mechanisms that could raise the inclinations 
of an initiall y low-i popula tion to values as large as 30° which have been observed in the 
Ku iper belt. iGomed (120031) investigated how the outward migration of Neptune proposed 
by iMalhotra ( Il993l . llQQsi ) could scatter objects from ~ 25 AU onto high-i orbits in what 
is now the classical Kuiper belt region. This study suggested that the high-z population 
formed closer to the sun and was emplaced in the classical Kuiper belt during planetary 
migration, whereas the low-eccentricity, low-inclination part of the classical Kuiper belt 
represents a primordial, relatively undisturbed population. It also led to the speculation that 
the observations of a correlation between in clinations and color in the classical Kuiper belt 
( jjewitt &: Luull200ll : iTrujillo &: Brownll2002l ) might be du e to different origins for the ob jects 
rather than environmental effects (see review article by iDoressoundiram et al.l (120081 )). In 
addition to the color-inclination correlation, there have been correlations p roposed between 
inclination and other phy s ical properties su ch as size, albedo, and binarity (ILevison fc Stern 
200 ll : iBrucker et ahlboOQl : INoII et al.l 120081 ). 



In this paper, we investigate the degree to which planetary perturbations may have 
caused dynamical changes in the inclinations of both the resonant and non-resonant KBOs in 
the classical Kuiper belt region over gigayear timescales. If the inclination-physical property 
correlations are real and are primordial, then the observed low- and high-i CKBOs cannot 
have undergone much dynamical change in inclination, else the correlations would be erased. 
If the differing physical properties are due to environmental effects (such as space weathering 
or collisions), the effect that changes in i could have on the correlations is less clear because 
it would depend on the timescale of the change in i and how quickly the environmental effect 
operates to change the properties of an object. In either case, the goal of this paper is to 
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quantify how much mixing there is hkely to be between the currently observed low- and high- 
i populations solely as a result of dynamical evolution under the current architecture of the 
solar system and to quantify how this might affect the observed correlations between physical 
properties and inclination. We outline our approach to modeling the current population in 
Section [2] and present the results of our numerical simulations in Section |3l In Section |4] we 
describe the dynamical mechanisms responsible for the changes in inclination for both the 
resonant and non-resonant populations, and we apply the results of our simulations to assess 
the correlations between physical properties and i. Section Ogives a summary of our results. 



2. Modeling the Classical and Resonant Populations 

To build a model of the classical Kuiper belt, we started with the observed set of KBOs 
with semimajor axes, a, in the range 33 < a < 50 AU and that have been observed at three 
or more oppositions, as these objects have more accurate orbit determinations than objects 
with fewer observations. Orbital elements and observational data for all objects used in this 
study were obtained from the Minor Planet Center (MPC) website0 on 2008 May 15; 598 
objects from this database met our above criteria. We generated 10 clones of each of these 
objects by varying the object's nominal determined orbit by small amounts (about one part 
in 10^, comparable to the uncertainties in the orbital elements) in a, eccentricity, e, argument 
of perihelion, oj, and mean anomaly. 

With the set of initial conditions obtained above, we performed a 10 Myr numerical 
integration of the orbital evolution under the current architecture of the solar system. The 
integration s were performed wi t h a m ixed variable symplectic integrator based on the al- 



gorithm of IWisdom &: HolmanI (jl99ll ). The KBOs were treated as massless test particles 
moving under the gravitational field of the Sun and the four outer planets; the mass of the 
Sun was augmented with the total mass of the terrestr ial planets. Initial c onditions for 



the planets were taken from the JPL Horizons service! ( Giorgini et al.lll996l ). We output 
the orbital elements every 5000 yr. We then checked for stability, for membership in mean 
motion resonances (MMRs) with Neptune, and we calculated the average values and rms 
variations of the orbital elements over this time period. Any KBOs for which the majority 
of the clones approach within a Hill radius of any planet during this 10 Myr period are 
discarded from our model population because such objects will generally leave the classical 
Kuiper belt population on short timescales after such an encounter; objects like these with 



^ http://www.cfa.harvard.Gdu/iau/mpc.html 
^http: / /ssd.jpLnasa.gov/?horizons 
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dyna mical lifetimes less than 10 Myr are best classified as Centaurs f lTiscareno fc Malhotra 
20031 ) and are therefore not part of the classical Kuiper belt. Only five of the objects were 
found to be unstable on Myr timescales. 

The membership in MMRs is determined as follows. An object is deemed resonant if its 
clones collectively spend more than 50% of the 10 Myr simulation with a librating resonance 
angle. While 50% is an arbitrary choice, in practice most objects are very obviously resonant 
(almost all of the clones librate for the entire 10 Myr simulation). The MMRs with Neptune 
(labeled a.s q : p MMRs) are characterized by resonant angles of the form 

4> = Q^kbo — P^N — fkbo'^kbo — f'N'^N — Skbo^kbo — Sn^N (1) 

where A, vu, and Q refer to the mean longitude, longitude of perihelion and longitude of 
ascending node, respectively, the subscripts kbo and refer to the elements of a KBO 
and to Neptune, respectively, and p,q,rkbo,rN, Skbo, sn are integers. Rotational invariance 
imposes the constraint q — p — Vkbo — — s^bo — Sn = 0. For objects with semimajor 
axes greater than Neptune's, the MMRs of interest have q > p > 0, and the order of 
the resonance is defined as \q — p\- For small e and i, perturbation the ory informs that the 
strength of an MMR is proportional to e|2°g''°'el^^l(sinzfc;,o)l'^''''°l(sinijv)''^^' ( Murray fc Dermott 
19991 ). Because rkbo + + Skbo + sn = q — P, the latter value is referred to as the order of 



a resonance, and it is a measure of the strength of the resonance. For small or moderate 
values of e and sinz, the lower order resonances are the stronger resonances. To simplify 
the process of checking for membership in the q : p MMRs, we limit ourselves to resonance 
angles of the form 

4> = Qhbo - P^N - (q- p)^kbo- (2) 



Previous works (j Chiang et al.l l2003l : lEUiot et al.l |2005| ) have found no examples of KBOs 
in an MMR where some other allowed combination of Vkbo^fN^ Skbo^ and s^r in Equation [T] 
showed libration and the above angle did not; a spot check of different resonance angles 
for our simulation supports this simplifying assumption. We identify 194 resonant objects 
including 108 objects in the 3:2 MMR, 28 in the 7:4 MMR, 16 in the 2:1 MMR, and 42 
amongst numerous other resonances; the highest order resonance identified is order in 
the KBO's eccentricity. The locations of the MMRs identified in this population, as well as 
the numbers of objects in each MMR, are shown in Figure [H Table [1] lists the designations 
of all the objects found i n resonances, noting objec ts n ot previously i denti fied as resonant 
bv lciadman et al.l J2008h . Ihvkawka fc Mukarj2007h . or lchiang et"al] J2003h . 



From the remaining non-resonant objects, we wish to build an orbital distribution that 
represents the intrinsic distribution for the classical belt by accounting for biases in the 
set of observed objects. Debiasing the observations requires knowledge of the discovery 
circumstances for each object, therefore in the non-resonant case we limit ourselves to objects 
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discovered by well characterized surveys; Table [2] g ives a complete lis t of su rveys used in this 
work. This criterion yields 278 CKBOs. Following I Vo Ik fc Malhotral fl2008[ ). the debiasing is 
achieved in our model by assigning each object's clones a weighting factor that is determined 
by the probability of detection for the object; KBOs with orbits that make them less probable 
to detect are given larger weights because observations thus far have a high probability of 
having missed many similar objects. KBOs with orbits that keep them nearly always within 
the observational range of telescopic surveys are weighted less because it is likely that almost 
all similar objects are within our observational data set. 

The primary assumption of our debiasing procedure is that an object's detectability 
only depends on the limiting magnitude and ecliptic latitude ranges of the observational sur- 
vey. This is a fair assumption for non-resonant objects but not for the resonant ones. The 
debiasing of resonant objects is more problematic because the probability of their detection 
depends not only on the limiting magnitude and ecliptic latitude of the observational sur- 
vey, but also on the longitude of the observations with respect to Neptune; the libration of 
the resonant argument restricts the object's longitude relative to N eptune and therefore the 
detection probability varies with longitude ( iKavelaars et al.l 120081 ). This is further compli- 
cated by the possibility of the Kozai resonance within the mean motion resonance (discussed 
further in section I4.1.2p . which restricts where the location of the KBO's perihelion will 
occur relative to its lo ngitude of ascending node, and therefore relative to the ecliptic plane 
( IKavelaars et al.ll2009l ). For example, Pluto is in both the 3:2 mean motion resonance and 
the Kozai resonance. The latter causes Pluto's argument of perihelion to librate around 
90 degrees. Physically, this means that the location of Pluto's perihelion librates around 
its maximum excursion above the ecliptic plane. This phenomenon adds to the bias against 
detecting such objects because they are furthest from the ecliptic plane when they are bright- 
est. Without knowledge of all the details of an object's resonant behavior, as well as a very 
specific pointing history for the observati onal survey, it is impossib le to correctly calculate 
a discovery probability (see, for example, IChiang &: Jordan! (120021 ) for a discussion of how 
discovery probabilities for the 3:2 and 2:1 MMRs depend on resonant behavior). This makes 
an accurate debiased model of their intrinsic orbital distribution beyond the scope of this 
paper; instead we separate the analysis of the resonant objects from the non-resonant and 
only apply our debiasing procedure to the non-resonant objects. 

To test our debiasing procedure, we compare our debiased inclination distribution to 
the debia sed inclination distr ibution reported for CKBOs from the Deep Ecliptic Survey 
(DES) by iGulbis et al.l ( l2010l ). Because the DES classification scheme differs significantly 
from our own, we only compare the sub set of our CKBOs that overlaps with their definition. 
Under the DES classification scheme (lElliot et al.ll2005l ). CKBOs are non-resonant KBOs 
with eccentricities < 0.2 and Tisserand parameters with respect to Neptune > 3. Using this 
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definition, iGulbis et al.l ( 120101 ) include 150 objects in their debiased CKBO i distribution that 
are also included within our CKBO dataset; the debiased inchnation distribution for these 



150 CKBOs is shown in Figure [2J Following iGuIbis et al.l ( 120101 ) . we smooth our debiased 
z-distribution at higher inclinations (where there are fewer observations on which to base 
the debiased distribution) by using variable bin sizes: 1° for i < 4°, 2° for 4 < z < 10°, and 
5° for i > 10°. The entire distribution is normalized to 1 and reported in Figure [2] as the 
fraction of CKBOs per degree in inclination. We performed a Kolmogorov-Smirnov (K-S) 
test to quantify the similarity of the two debiased inclination distributions. The K-S test 
measures the maximum difference, D, between the cumulative distributions of two individual 
data sets an d then evaluates t he probability that the two are drawn from the same parent 
distri bution (iPress et al. ll992[ ). The K-S test comparing the debiased CKBO i distribution 
from iGulbis et al.l (120101 ) to our debiased i distribution for the 150 CKBOs in our sample 
that overlap with the DES sample yields a 65% probability that the two distributions are the 
same. Compared to the DES debiased distribution, our distribution is slightly systematically 
skewed toward lower inclinati ons. This difference is due to our simplifying assumptions in the 
debiasing procedure: whereas iGulbis et al.l ( l2010l ) use the complete details of each discovery 
image (such as limiting magnitude and ecli ptic latitude) to ca lculate debiasing factors, we 
used the DES survey averages described bv lEUiot et al. J2005I). Given this limitation, our 
debaised distribution is close enough to that of ICulbis et al. (2010) to give us confidence in 
our debiasing procedure. A further test of our debiasing procedure is to compare th e resulting 
inclination distribution to that found for the classical Kuiper belt by iBrownl (j200l[ ). who used 
a definition of CKBOs similar to our own; we discuss this comparison in Section [31 

The set of 194 resonant and 278 non-resonant objects in the classical belt region selected 
from the 10 My simulation were then integrated for 4 Gyr with output every 1 Myr to study 
the long-term orbital evolution of each object under the perturbations of the four outer 
planets. As in the case of the 10 Myr integrations, 10 clones of each object were integrated 
as massless test particles, and we terminated the integration of any particle that approached 
within one Hill sphere of any planet during the 4 Gyr simulation. In the analysis in the 
following section, we refer to the particles that survive to the end of the 4 Gyr simulation 
as "stable", and the others (that have dynamical lifetimes less than 4 Gyr) as "unstable". 
The vast majority of the unstable particles in our sample have gigayear lifetimes, and are 
thus expected to be representative of the presently existing Kuiper belt; we include their 
contributions in our time-averaged results. 
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3. Analysis and Results 



To analyze the inclination evolution of our test particles, we must choose an appro- 
priate coordinate system for our results. The standard output for the numerical integra- 
tions are heliocentric orbital elements; in our case these are referenced to the J2000 coor- 
dinate system, meaning all inclinations are measured relative to the ecliptic plane of the 
epoch. Because we are interested in changes in the inclinations of the test particles, it is 
preferable to measure those inclinations relative to a more physically meaningful reference 
plane; the best reference plane for determining changes in the orbits of the KBOs would 
be the mean plane of the Kuiper belt. Several works have calculated t he mean plane of 
the Kuiper be . 



t based on observations an d secular perturbation theory (lElliot et al.l 12005 
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2008|); this is 



sl); this is not a desirable feature for a reference plane for 



our numerical simulations. A more conveni ent and easily c a lculat ed re ference plane is the 
invariable plane of the solar system. Both IChiang fc Choil ( 120081 ) and lEUiot et al.l (j2005[ ) 
find only smal l devia tions between the invariable plane and the average Kuiper belt plane; 
Brown fc Paril (120041 ) found a larger difference between these planes, but the difference was 
still less than 1°. We will therefore reference all KBO inclinations in this paper to the 
invariable plane. To do this, we calculate the location of the barycenter for the Sun and 
the four outer planets, and then we calculate the total angular momentum vector for Sun 
and outer planets relative to the barycenter; the plane normal to this vector becomes our 
reference plane for the orbital inclinations. Adopting the invariable plane diminishes the 
amplitude of the short-term oscillations in the KBO test particles' inclinations compared 
to when the ecliptic plane is used. This is because their short-term inclination evolution is 
driven by the short-term inclination evolution of the planets. If the ecliptic plane is used 
as the reference plane, the planets' inclinations average 1.5 — 2° over Myr timescales; if the 
invariable plane is used instead, this average is 0.3 — 1°. The reduction in the amplitude of 
the corresponding changes in the test particles' inclinations is illustrated in Figure [3l which 
plots the rms variation of each test particle's inclination vs. the mean inclination for the 10 
Myr simulation. This reduction in the 'typical' spread in the inclinations expected over Myr 
timescales makes it easier to distinguish the test particles with larger inclination changes in 
the longer simulation. Unless noted otherwise, all inclinations in this paper are referenced 
to the invariable plane. 

We can use several measures to determine the mobility of the KBOs in inclination space. 
Figures m and [5] show the rms variation in inclination over the initial 10 Myr simulation and 
the 4 Gyr simulation, respectively, for all the test particles as a function of average semimajor 
axis over the same time period (note the different scales on the y-axes of Figures H] and [5]) . 
The stable test particles (those that survive the entire 4 Gyr simulation) are shown separately 
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from the particles that have close encounters with Neptune. The inclination variations in 
the 4 Gyr simulation are generally larger than in the 10 Myr simulation, especially for the 
resonant objects. The maximum changes in any one particle's inclination over the entire 
simulation can be many times these rms values. Figure E] shows the maximum Ai for 
objects that survive the entire simulation and those that do not, where both groups are 
separated into CKBOs and resonant objects. Not surprisingly, the unstable objects undergo 
the largest changes in i, but several of the stable CKBOs (less than 1%) and quite a few 
of the stable resonant objects (8%) show changes of more than 10°. Of the total debiased 
CKBO population, 10% of objects deviate from their initial 10 Myr average i by more than 
5° at some point during the simulation (2% of the stable CKBOs and 84% of the unstable 
CKBOs), and 4% deviate by more than 10° (<1% of the stable CKBOs and 40% of the 
unstable CKBOs). For the observed population of resonant objects, these percentages are 
51% (38% of the stable objects and 87% of the unstable objects) and 23% (8% of the stable 
objects and 64% of the unstable objects), respectively. 

Given that some objects' inclinations are quite changeable over time due to planetary 
perturbations, it is important to examine the stability of the overall inclination distribution. 
Figure [7] shows snapshots of the inclination distribution for the debiased non-resonant CK- 
BOs at the beginning of the simulation, 500 Myr into the simulation, and at 4 Gyr. (Note 
that we use the same binning in this figure as in Fig. 2.) The results of a K-S test (described 
in Section [2]) show that there is a 52% probability that the initial and final snapshots are 
drawn from the same inclination distribution, and there is a greater than 90% probability 
of similarity between all snapshots taken after 500 Myr. The lower K-S probability for the 
early snapshots compared with the final snapshot is most likely due to the circumstance that 
if an object on an unstable orbit is heavily weighted in the distribution and is then either 
removed from the simulation or evolves onto a more stable orbit, the resulting change in the 
overall distribution can be statistically large. The difference between the initial and final 
distribution is not large, and a K-S probability of ~ 50% is not so small as to rule out 
the possibility that the distributions are the same; the 90% K-S probability after the first 
500 Myr shows that, even if some objects undergo large changes in inclination, the overall 
distribution is quite stationary. 

The debiased inclination distribution we find for the CKBOs (Figure [7]) shows a distinct 
peak at low-inclinations {i < 5°) with a substantial part of the population spread over a wide 
range of higher inclinations. The rea der will note that th is is significantly different from the 



i-distribution for the DES CKBOs by lGulbis et al.l (120101 ) (see Figure [2]): the notable paucity 



of high-z CKBOs in Figure [2] compared to Figure [7| is owed entirely to the differing definitions 
for CKBOs. The cutoff in the Tisserand parameter with respect to Neptune (T > 3) in the 
DES definition is the major contributor to the difference; for example, a KBO with semimajor 
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axis a = 40 AU and an eccentricity of 0.1 must have an inclination less than ~ 12° and a 
KBO with a = 45 AU and an eccentricity of 0.1 must have an inclination less than ~ 16°, 
to qualify as a CKBO using the DES definition. Th e eccen t ricity cutoff (e < 0.2) adopted 
by DES also excludes a few of our high-z CKBOS. iBrowru (120011 ) use a CKBO definition 
much more similar to ours than that of the DES, the only differences being the adoption of 
a slightly narrower range in semimajor axes (40 < a < 48 AU) and the exclusion of a few of 
the higher eccentricity KBOs in that semimajor axis r ange for havin g perihelion distances 
very similar to that of scattered disk objects. Following iBrownl (120011 ). we model our initial 
inclination distribution to a sum of two Gaussians multiplied by sin i 



smz 



Aexp 



+ (1 - A) exp 



(3) 



We find the following best-fit parameters: A = 0.95±0.02, ai = 1.4±0.3°, and a2 = 15±3°. 
The se may be com pared with A = 0.93 ± 0.03, ai = 2.2lQg, and (J2 = 17 ± 3° from Table 
1 of iBrownl (120011 ). The largest difference between the two best-fits is in the width of the 
narrow component; this difference is mainly owed to the differing reference planes; iBrown 
(120011 ) used the ecliptic plane as the reference plane for the inclinations, while we use the 
invariable plane. Our fit indicates that 83% of the CKBOs ar e in the wide G aussian and 17% 
in the narrow Gaussian, similar to the 81% and 19% found by IBrownl (120011 ). The inclination 
at which the narrow and w i de com ponents have equal numbers of objects is i = 5° in our 
model and i = 7° in IBrownl (120011 ). Both fits are plotted as the smooth curves in Figure [71 
While the data appear to support the existence of both a \ow-i and a high-z component in 
the inclination distribution. Figure [7] shows that the model curves do not fit the data very 
well. This poor fit is at least partially due to the noisy nature of the high-z portion of our 
debiased distribution: there are still relatively few observed objects at large-i. The poor fit 
could also be an indication that a double Gaussian is not a very good model for the intrinsic 
distribution. Gulbis et al. (j2010 ) try various other functional forms for fitting the debiased 
DES inclination distribution and find that a Gaussian plus a generalized Lorentzian provides 
an improved fit; how ever, the impr ovement over the double Gaussian is small, so for ease of 
comparison with the IBrownl (120011 ) results we only report a fit to the double Gaussian form 
of the inclination distribution. 

From the 4 Gyr simulation we want to quantify how well an object's current inclination 
correlates with its past or future inclination; i.e. what is the probability that a currently 
'low' inclination object has always been in the 'low' inclination group? The answer to this 
question may depend on how 'low' inclination is defined. We adopted a cutoff inclination 
value, and then examined how the degree of mixing between the high- and low-inclination 
populations depends on the value of the cutoff inclination. Figure [8] shows the time-averaged 
fraction of all non-resonant and resonant CKBOs that can be found at inclinations that are 



- 10 - 



inconsistent with their initial (i.e. current epoch, 10 Myr average inchnation) classification 
into either the high- or low-z populations as a function of the i cutoff that defines the low- 
and high-inclination populations. For example, fo r a cutoff of 5° (which is soraetimes used to 



distinguish the so-called 'hot' and 'cold' CKBOs ( iDoressoundiram et al.ll2002t iBrucker et al. 



20091 )) one would expect 1.5% of all non- resonant CKBOs and 11% of resonant objects 
to currently be at inclinations inconsistent with their primordial inclinations, because the 
populations as a whole spend 1.5% and 11% of their time at inclinations inconsistent with 
their current classification^. Our results in Figure [H] also show that the inclination of ~ 5° 
is justified as a divide between the hot and cold population as it corresponds to a local 
minimum in the 'misclassification' distribution. 

We note that if the ecliptic plane were used for this calculation instead of the invariable 
plane, the result is more inclination mixing (and therefore a fuzzier distinction) between the 
two classes, as shown in the second panel of Figure [H For both the CKBOs and the resonant 
KBOs, the percentage of objects we expect to find at inclinations inconsistent with their 
initial classification is smaller than the percentage of objects found to have deviated from 
their initial inclinations by large amounts at some point in our simulation. This is because 
objects that experience large deviations in i do not spend a large fraction time at the extreme 
values. So even though half of the resonant KBOs will at some point have an excursion more 
than 5° away from their initial i, only 12% of these objects spend more than 10% of their 
time more than 5° away from their initial i. Similarly, only 6% of CKBOs spend more than 
10% of their time at values of i more than 5° away from their initial i. 



4. Discussion 
4.1. What drives the changes in inclination? 

4-1.1. CKBOs (non-resonant) 

Approximately 10% of the total debiased CKBO population experiences a change in i 
of more than 5° during our 4 Gyr simulation. Here we examine the causes of these relatively 
large changes in inclination to better understand how the inclinations of CKBOs evolve over 
time. To do this we examined the time histories of the test particles exhibiting the largest 



■^If. instead of using the average inclination from the 10 Myr simulation, we use the instantaneous incli- 
nation with respect to the invariable plane, the percentage of misclassified objects using a 5° cutoff between 
low- and high-i populations is still 1.5% for the CKBOs, but for the resonant objects the percentage increases 
to 13%. 
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overall changes in inclination as well as those exhibiting small changes. Objects with stable 
orbits and inclinations tend to undergo variations in z of ~ 1° that occur on timescales 
comparable to or shorter than our output resolution of 1 Myr (an example of this is shown 
in Figure [H]). This is consistent with what is expected from the secular perturbations of the 
planets on a massless test particle; near Neptune, the plane of a KBO's orbit will wobble 
about the invariable plane with an amplitude s imilar to Neptune's proper inclination; th e 
latter is ~ 0.7° inclined to the invariable plane (IBrown fc Panll2004l : IChiang fc Choill2008l ). 

Most (about 80%) of the classical belt test particles with large Ai have changes in i 
that can be attributed to distant encounters with Neptune. The time resolution of our 4 
Gyr simulation output is not fine enough to detect the encounters, but we can diagnose 
such encounters with Neptune because they manifest as discrete jumps in the semimajor 
axis, eccentricity, and inclination of the test particles, and these occur when the particle is 
at a local minimum in the time history of the perihelion distance. An example is shown in 
Figure [TOl where an ~ 8° increase in i occurs at a perihelion distance q of about 32 AU. 
These moderately distant encounters with Neptune certainly explain the large changes in i 
for the unstable subset of CKBOs. The values for the typical spread in i and the maximum 
Ai reported here exclude the 50 Myr period prior to the closest approach with Neptune 
(within one Hill radius) that removes a test particle from our simulation, so the variations 
in i shown in Figures [5] and [6] represent changes that are occurring due to more distant 
encounters. Unstable test particles can wander substantially in inclination over timescales 
much longer than 50 Myr before having an encounter close enough to completely remove 
the particle from the classical belt region. Approximately one-third of the test particles 
experiencing large Ai as a result of a distant encounter with Neptune persist in the classical 
belt for timescales longer than 250 Myr after the encounter. A small subset of the stable 
CKBOs will also experience distant encounters with Neptune that lead to significant changes 
in i without leading to a close enough encounter to completely remove the test particle (at 
least not within the 4 Gyr of our simulation). 

Approximately 20% of the test particles (both stable and unstable) with large maximum 
Ai values in the classical belt owe their inclination variability to brief stints of libration within 
an MMR. Careful examination of Figures O and [H] shows that some of the CKBO test particles 
in close proximity to resonant particles (as measured by their average semimajor axis values) 
also have large rms variations in i and/or maximum Ai, similar in magnitude to the nearby 
resonant particles. Some test particles are captured for long enough in MMRs that we detect 
libration of the resonance angle in our simulation output, and we have identified temporary 
capture in most of the occupied resonances and a few very high-order MMRs with Neptune 
such as the 19:10 and 21:11. Short-lived resonance libration is difficult to ascertain for test 
particles that spend less than a few tens of millions of years in resonance, due to the lower 



-VI- 



time resolution of our simulation out put. Temporary capture i nto various resonances in the 



classical belt region was also noted by lLykawka fc Mukail (j2006[ ) in a long-term simulation of 



test particles started near the 7:4 MMR. In our simulation, we find that some test particles 
start in a resonance, but exit the resonance without being removed from the classical belt 
population. Some of the resonant test particles also experience only intermittent libration 
(see Section 14.1. 2p . so it is possible that objects we classify as non- resonant are close enough 
to a resonance boundary to experience occasional librations. Such chaotic diffusion across 
resonance boundaries is the cause of many of the cases of large inclination variability of 
CKBOs. 

A few percent of the stable particles with large in the 4 Gyr simulation cannot 
be explained with either of the above two mechanisms - distant Neptune encounters or 
temporary resonance capture. There is a secular resonance at ^ 41 AU that ha s been 



shown to raise the inclinations of initially low-i test particles (IKnezevic et al.lll99ll ). and, 
as discussed above, there are also many currently occupied MMRs in the 40 to 45 AU 
range which might be expected to influence the orbital evolution of nearby test particles. 
Because these dynamical mechanisms work mainly in the inner part of the classical belt 
region (a < 45 AU), we expect the amplitude of changes in i to decrease in the outer region, 
with the exception of test particles near the 2:1 MMR. In the 10 Myr simulation, this is 
indeed what we find; the typical spread of i for the stable objects (as seen in the left panel 
of Figure Hj) is small in the outer region, and it increases moving inward toward both the 
5:3 MMR at ~ 42.3 AU and the secular resonances at 41-42 AU. However, in the longer 
simulation (Figs. E] and E]) we see some particles with surprisingly large rms variations in i 
and maximum changes of 5 — 10° in the semimajor axis range of 45 to 47 AU where there 
are no strong MMRs with Neptune and no known secular resonances. This zone does have 
several higher order mean motion resonances, such as Neptu ne's 11:6 and 7:19, as w ell as 



a three-body resonance involving both Neptune and Uranus ( iNesvorny h, Roig]|200ll ). We 
looked for librations of candidate resonance angles; our output cadence of 1 Myr (for the 
4 Gyr simulation) can only detect libration behavior on timescales greater than ~ 10 Myr; 
we did not detect any such long-lived librations. For our test particles between 45 and 47 
AU, the effects of chaos arising from the overlap of these high-order interactions may be 
the cause of the larger variability in i that we find in the 4 G yr simulation (bu t do n ot 



find in the 10 Myr simulation). This conjecture is supported by iNesvorny &: Roigi (120011 ) 's 
finding of a small positive Lyapunov exponent in this range of semimajor axis, although their 
analysis was confined to very low inclination, low eccentricity orbits. A deeper analysis of 
the dynamics of objects in this zone is left for a future investigation. 



There still remain a few test particles notable for their peculiar longer-timescale large 
inclination variations for which we have not identified any causative dynamical mechanism. 
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A typical example of these is shown in Figure [TT] i varies from 2° to 9° with a periodicity 
of about half-a-gigayear (with a corresponding synchronized variation in eccentricity from 
0.06 to 0.15). Test particle clones of four other CKBOs show similar orbital element time 
histories, but with shorter periodicities ranging from about 50 Myr to about 250 Myr. These 
timesc ales are much longer than the libration time scales of MMRs which are typically 10^- 
10^ yr (jMalhotralll996l : iTiscareno fc Malhotrall2009l). or the time scale of the Kozai resonance 
within MMRs which is typically 10^-10^ yr f IChiang et al.ll2003f l. The semimajor axis of the 
test particle is also inconsistent with the location of t he fgj^^is and secular resonances 
which are all located at 40-42 AU ( iKnezevic et al.lll99l[ ). We performed additional numerical 
simulations with the test particles exhibiting this peculiar behavior and were able to rule out 
numerical artifacts (we varied the integration timestep and integrated in both the ecliptic 
and invariable coordinate frames) as a potential cause. We also tested for membership in 
all of the nearby r esonances with Neptune, Uranus, and combinations of both planets that 
were identified by iNesvorny fc Roigl (j200l[ ). but found no librating resonance angles. We 
speculate that the cause may be some other very high orde r MMR we did not c onsid er, or 
possibly a super-re sonance such as identified for Pluto by IWilliams fc BensonI ( 1l97ll ) and 



Milani et al.l (119891 ): verification of this is left for future work. 



4.1.2. Resonant KBOs 

The resonant objects in our simulation generally show more variation in i than the 
non-resonant objects do; large changes in i (Ai > 5°) occur for 51% of the resonant objects 
over 4 Gyr compared to only 10% of non-resonant objects (see Section |3]). The variations 
in i for the resonant KBOs occur on multiple timescales, with larger changes occurring on 
longer timescales; this can be seen by comparing Figures H] and O Libration within an MMR 
causes variations in the orbital elements over timescales similar to the libration period of 
the resonance, which is typically a few tens of kiloyears for the resonances in this study. 
Resonant particles also experience the same secular variations in i on megayear-timescales 
that we describe above for the CKBOs. Over longer timescales (a few tens of megayears) , the 



Kozai resonance (jKozailll962l ) accounts for many of the resonant particles having maximum 
inclination excursions > 5°. The Kozai resonance is characterized by the libration of the 
argument of pericenter, and the preservation of the Kozai integral, 

e = V(l -e2)cosi (4) 

The effect of this resonance is large anti-correlated oscillations in eccentricity and inclination. 
We find that many test particles do not librate in the Kozai resonance for the entire simu- 
lation, but instead experience intermittent libration which results in coupled large changes 
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in e and i. We also find that the resonance angle (j) (as defined by equation [2]) can exhibit 
intermittent changes in amplitude (correlated with intermittent Kozai libration) which also 
affects the inclination evolution. A typical example of this type of evolution is shown in 
Figure [T21 This type of intermittent libration of both cj) and the argument of perihelion 
has been fo und in previous work ex ploring the dynamics of the 7:4 MMR in the classical 



belt region (ILykawka fc Mukaill2005l ). This behavior is responsible for the stable resonant 



population having much larger maximum values of Ai than the stable non-resonant CKBOs. 

We also note that some members of the 7:4 and 9:5 MMRs diffuse out of resonance and 
evolve into stable orbits in the classical belt region. (These objects are indicated in Tabled]). 
This phenomenon likely contributes to the population of CKBOs that exhibit larger changes 
in i (as discussed in Section [4. l.ip . Objects like these 'resonance escapees', whose dynamics 
at early times were dominated by the effects of MMRs, but later cease to librate in the 
resonance, are probably present in the observed sample of currently non-resonant objects. 



4.2. Reassessing the correlations between inclination and physical properties 



Here we examine the proposed correlations between physical properties and inclination 
using the results from our numerical simulations to represent the time-varia bility of the 



inclinations of the observed objects. We use the Spearman rank correlation test (jPress et al. 



19921 ) as a quantitative measure of the correlations between z and each physical property. This 
test has been used in previous studies of this problem (see iDoressoundiram et al.l (120081 ) and 
references therein); it is chosen because the results of the test do not depend on any assumed 
form for the correlation (i.e. the relationship need not be linear). Values of Rg close to ±1 
indicate a strong correlation or anti-correlation between the two parameters tested. The 
significance level of the correlation is reported as the number of standard deviations above 
or below zero (no correlation) the value of Rg lies compared to the null hypothesis that the 
two variables are uncorrelated. For each physical property, we perform the correlation test 
using the average value of i for the objects from our 10 Myr simulations, which is a better 
representation of their proper inclinations than the instantaneous observed value of i; we use 
the rms Ai from the 10 Myr simulation as a measure of the time variability of i (using the 
average i and rms Ai values from the 4 Gyr simulation does not alter the conclusions that we 
describe below). The physical properties we consider for the correlation test are the spectral 
gradient and absolute magnitude; we do no t consider albedo because only ~ 30 CKBOs 
have measured albedos (IBrucker et al.ll2009l ). We compare inclination distribution of the 
binary CKBOs to that of the apparently si ngle KBOs, as well as the i distribution of large 
CKBOs to the distribution determined by iBrownl (120011 ) for the classical belt as a whole. 
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These comparisons are done using the Kolmogorov-Smirnov test (described in Section [2]) to 
measure the probabihty that the inchnation distributions being compared are the same. 



4-2.1. Spectral Gradient 

The colors of KBOs are assigned by photometric characterization of their reflectance 
spectra in various broadband filters (BVRI for the visible wavelengths). The commonly used 
color indices (B-V, V-R, etc) represent the differences in the observed magnitude of the KBO 
in the two filters. An additional color parameter that is determined from these broadband 
filter measurements is the spectral gradient, S; this is a measurement of the reddening of 
the KBO's reflectance spectrum over the wavelength range considered, usually expressed as 
percent reddening per 100 nm. We use the spectral gradient in our color-z correlation tests 
because S can be calculated for any KBO that has been observed in at least two visible 
broadband filters, and it has been shown to co rrelate with the commonly measured V-R 



color for KBOs (see lDoressoundiram et al.l (120081 ) for a full discussion of this). The results of 



our statistical tests do not depend on this choice (results for B-V, V-R, and B-R color indices 
are very similar to those for S), but using the spectral gradient provides us with the largest 
possible observational data set. Figure [13] shows the spectral gradient vs. 10 Myr average 
inclinations for 80 CKBOs. All values for S are taken from the MBOSS d atabas jll maintained 



by O. Hainaut of the European Southern Observatory and described by iHainaut fc Delsanti 



(I2OO2I). 



In our calculations of the correlation coefficient, Rs, using the Spearman rank test, we 
want to include information about the observational uncertainties for the spectral gradient 
and the time variability of the inclinations (as represented by the rms variations in i cal- 
culated from the numerical simulations). Including the uncertainties for both the physical 
property and the inclination is achieved by randomly assigning values to both the observed 
physical property and i, assuming that each variable has a normal distribution; this is re- 
peated 1000 times and Rg is calculated for each permutation of the data set. The resulting 
range and significance of the values of Rg yields a measure of how the observational errors 
in S and the dynamical evolution of i affect the correlation. Using the average i and its rms 
value from the 10 Myr simulations as well as the observational errors in S, this procedure 
yields Rg = —0.46 ±0.07 (4±0.5cr significance). This can be compared to Rg = —0.52 (4.6cr 
significance) for the test when we use only the best-estimate values of S and the average 
values of i. Most of the reduction of Rg and its associated significance is due to the observa- 



"'http: / / www.eso.org/ ^ohainaut /MBOSS / 
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tional errors in S. When we account for the uncertainties in S but do not take account of the 
time-variabihty of i, we find Rg = —0.47 ± 0.05 (4.2 ± 0.5a significance). This is consistent 
with the fact that most of our observed sample of CKBOs with measured values of S are 
stable and have fairly low variations in inclination. 

Another factor that could affect the correlation between spectral gradient a nd inclination 
i s the existence of the Haumea (2003 ELgi) collisional family discovered by iBrown et al. 



( 120071 ). Haumea family members share a distinctive water ice absorption feature, and they 
are clustered in orbital element space with inclinations near ~ 25°; the family m embers 



also tend to have neutral colors ( ISchaller fc Brown 



2008 



Ragozzine Sz Brownll2007D. There 



are six Haumea family members (as defined by iBrown et al.l (l2007f ) and iRagozzine fc Brown 
( 120071 )) in our set of CKBOs that have measured spectral gradients: Haumea (2003 ELgi), 
1995 SM55, 1996 TOee, 1999 OY3, 2002 TX300, and 2003 UZny. These six objects comprise 
one third of the CKBOs with i > 25° that have observed spectral gradients, and all six 
have neutral spectral gradients, creating a cluster in the plot of S vs. i (Figure [T3|) . To 
test how this cluster of related objects affects the correlation coefficient for S and i, we 
calculated Rg for the CKBOs without these six objects in the data set. Removing the 
Haumea family members, but not taking the uncertainty in S, or the time- variability of i 
into account, we find Rs = —0.42 with a 3.6cr significance (compared to Rg = —0.52 with 
a 4.6(j significance with the Haumea family included). If we remove the family members 
and account for the time-variability of i and uncertainties in S, we find Rg = —0.36 ± 0.07 
(3.1 ± 0.6cr significance). We conclude that the neutrally colored, high-inclination Haumea 
family is partially responsible for the reported correlation between color and inclination in 
previous studies. 

The degree of correlation we find between color and inclination for the classical objects 
is smaller and less statistically significant than previou sly reported iri the l iterature. Several 
previous studies have found Rg between —0.6 and — 0.7 ( Jewitt fc LuiJlioOl : Doressoundiram et al, 



2002 



Trujillo fc Brownll2002t iPeixinho et al.ll2008l ): the most recent of these (jPeixinho et al. 



20081 ) analyzed the colors of 69 CKBOs, and reported Rg = —0.7 with a greater than 8a 
significance. In our analysis, even without accounting for the time- variability of i or the 
observational errors in S, the inclusion of 11 newly observed CKBOs in the analysis in this 
work (most of which have large inclinations) lowers the value and significance of the corre- 
lation coefficient to Rg = —0.52 with 4.6cr significance. This is a strong indication that the 
previously reported, Rg = —0.7 with 8a significance, did not reflect the true correlation value 
and its significance for the intrinsic population; if the significance level of the correlation were 
truly that high, it would have been extremely unlikely for additional observations to change 
the measured value of Rg by such a large amount. The large change in the significance of the 
correlation is likely a result of the observational incompleteness of the available sample of 
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CKBOs, especially at higher inclinations. Likewise, the inclusion or exclusion of the Haumea 
family members has a large effect on the significance of the correlation {Aa significance with 
the family, 3a without). Therefore, while the trend of less spectral reddening with increasing 
inclination appears to be significant at the 3a level in the CKBO dataset, this result should 
be interpreted with caution until the addition of new observations does not drastically alter 
the results of the Spearman rank test. 

We also tested the resonant KBOs for correlations between S and i, but we found no 
statistically significant correlation {Rg = —0.2 ±0.1 sig nificant at a < la level), whic h is con- 



sistent with previous work on the color-z relationship fiDoressoundiram et al.ll2008l ). Unlike 
the case for the CKBOs, for the resonant objects the errors in S contribute equally with the 
time- variability of i to decrease both the value and significance of Rs when compared to the 
use of only the best-estimate values of S and i. This is consistent with our numerical sim- 
ulations, in which resonant objects undergo much larger changes in i than CKBOs. Even if 
there were initial differences in spectral color between the low-z and high-i resonant popula- 
tions, subsequent dynamical evolution would tend to erase the trend through the significant 
inclination variability of resonant KBOs that we found in our dynamical simulations. 



4-2.2. Absolute Magnitude 



Here we examine the relationship between inclination and absolute magnitude. (In the 
absence of albedo measurements for most of the KBOs, their sizes remain uncertain, so 
absolute magnitude is used as a proxy for size). The question of how these two properties 
might be related can be framed in several different ways: (1) Do the low-i and high-z pop- 
ulations have different size distributions? (2) Do the large and small KBOs have different 
inclination distributions? (3) What is the degree of correlation between the two parameters 
(the Spearma n rank test)? These are different questions, the answers to which have different 
implications. iBernstein et al.l ( 120041 ) addressed the first question and detected a difference 
between the size distributions of the \ow-i and high-z KBO populations; the implications 
this has for the accreti onal histories of the se two populations has been explored in several 
works, most recently by iFraser et al.l (j2010[ ). In the present work, we focus on the latter two 
questions; the differences in the inclination distributions of the population of large and small 
KBOs might tell us if the two groups experienced different dynamical histories. 

There are some indications that the largest objects in the classi cal belt have an inclina- 
tion distribution that is different from that of the smaller objects; iLevison fc SternI (j200l[ ) 
report a statistically significant difference (97% confiden ce level) betweeii objec ts with ab- 
solute magnitudes H < 6.5 and those with H > 6.5. iLevison fc Stern! ( 1200 ll ) dismissed 



- 18 - 



observational biases as the source of the difference in inchnation distributions for large and 
small objects because the set of observed objects at that time did not show any correlations 
between size and ecliptic latitude at the time of discovery; discoveries of new KBOs at that 
time were dominated by ecliptic surveys. This is no longer the case as off-ecliptic surveys 
now account for many KBO discoveries, especially at brighter magnitudes. The observational 
surveys most capable of detecting faint objects have been performed near the ecliptic where 
there is a strong preference for finding low-inclination objects; the wider latitude surveys 
have brighter limiting magnitudes, and are therefore not as likely to detect small objects. 
The result is that observations of the brightest objects have fairly evenly sampled the in- 
clination distribution of the brighter objects, while the observations of the fainter objects 
are much more heavily sampled at low inclinations. Any test we use to compare the two 
distributions (one complete, and one biased) will show a difference between the distributions 
because of this sampling bias. A better way to test whether the largest objects have the 
same inclination distribution as the smaller objects is to compare the i distribution of the 
largest CKBOs to a debiased i distribution for all CKBOs. Here we will make the comparison 
using both the best-fit model to our debiased inclination distribution ( the init i al dist ribution 



shown in Figure [7]) and the best fit model distribution calculated by iBrownl (120011 ) for the 
classical belt. 

To minimize the effect of observational biases on the correlation analysis, we define 
the group of large CKBOs by adopting an absolute magnitude cutoff that yields a set of 
objects that is observationally fairly complete. Wide area surveys have detected most ob 



jects with apparent m agnitudes R < 21 within 30° of the ecliptic (ITrujillo fc Brown! 12003 



Schwamb et al.l 120091 ). The apparent magnitude R ^ 21 corresponds to an absolute magni- 
tude of ~ 5 at 40 AU, so we use this as the cutoff for observational completeness. At this 
absolute magnitude there is a sharp drop off in the fraction of objects detected at large eclip- 
tic latitudes: 60% of objects with H < 5 were detected at ecliptic latitud es between 5 and 



30° w hile only 15% of objects with 5 < H < 6 found at similar latitudes. iLevison fc Stern 



(120011 ) reported that CKBO's with H < 6.5 tend to have higher inclinations; if that trend is 
real, then the conclusions should hold for our smaller H cutoff (larger minimum size) . There 
are 26 CKBOs with H < 5 in our CKBO sample (even with a more stric t size cutoff, this 



i s a si gnificant improvement over the 8 objects with H < 6.5 available to iLevison fc Stern 



(120011 ) for their analysis). In Section 3, we noted that a natural cutoff inclination between 
the low and high inclination populations is z = 5° for our best-fit double gaussian model, and 
2 = 7° for Brown's 2001 model. Of the 26 large CKBOS, only one has i < 5° and three have 
i < 7°. We can compare these observed numbers of objects to the expected number of low-z 
CKBOs. Using our best fit parameters and their la uncertainties, we expect between 4 and 
8 objects with z < 5° in a set of 26 observations. The probability of observing only one such 
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object is 4_3 g%. Using the iBrowru (120011 ) best fit parameters and la uncertainties, we would 
expect between 5 and 8 CKBOs with i < 7° compared to the 3 we observe. The probabihty 
of this happening by chance is 1^1%. As noted in Sectional neither of these two model fits 
is a very good description of our debiased inclination distribution, but they do offer some 
description of the expected ratios of high- to low-inclination CKBOs. The subset of CKBOs 
with < 5 contains fewer low-i objects than expected if the true inclination distribution for 
this brightness range contains both high- and low-i Gaussian components; from the above 
analysis, we can estimate that there is a 4 — 7% probability (0.4 — 19% when we account for 
the uncertainties in the model fits) that this observation arises by chance. Considering our 



full ra nge of uncertainties, this result is consistent with the 3% reported by lLevison fc Stern 



fl200lh 



If the inclination distribution varies with size, we also might expect to see a correlation 
between H and i. For the smaller CKBOs, with i7 > 5, the correlation coefficient between H 
and average i from the 10 Myr simulations we find is Rs = —0.07 (1.3cr significance); when 
we include the effect of inclination variability from the 10 Myr simulations (as described 
in the previous section), we find Rg = 0.07 ± 0.01 (1.3 ± 0.3cr significance). This result is 
consistent with no correlation between H and i for H > 5. For the entire set of CKBOs, the 
correlation coefficient drops to Rg = 0.04 ±0.09, which is also consistent with no correlation. 

For the larger CKBOs, with H < 5, the correlation coefficient between H and average i 
from the 10 Myr simulations is Rs = —0.3 {1.5a significance); including the time variability 
of i, we find Rg = —0.29 ± 0.02 (1.5 ± O.la significance). These relatively low values of \Rs\, 
together with the low ~ 1.5a significance level of the Spearman rank correlation means that 
we do not have a definitive answer to whether H and i are correlated such that there is a 
different inclination distribution at different sizes of the KBOs. The anti-correlation between 
H and i is not highly significant for the large objects. There are fewer large objects than 
expected at small i, but there is a non-negligible ~ 4 — 7% probability that the discrepancy 
is a result of the relatively small number oi H < 5 objects rather than a true difference in 
the inclination distributions. Future wide-area surveys with deeper limiting magnitudes will 
allow us to compare the i distribution for a larger observationally complete set of bright ob- 
jects to the classical belt inclination distribution, at which time a more definitive judgement 
can be made. 

For the resonant KBOs, we find no significant correlation between H and i for the 194 
observed objects. Employing the same method for including the time-variability of i into the 
calculation of the correlation coefficient as in Section 14.2. H we find Rg = (—2 ± 2) x 10~^, 
consistent with the two parameters being unrelated. 
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4-2.3. Fraction of Binary Objects 



Noll et al.l ( 120081 ) searched 101 known CKBOs with inclinations ranging from 0.6 — 34° 
for binary companions and found a binary fraction of ~ 30% for objects at z < 5.5° com- 
pared to a binary fraction of ~ 10% at higher inclinations; if the comparison is restricted 
to the fraction of objects that are same-size binaries, the binary fractions are ~ 30% and 
~ 2% respectively. They used the Kolmogorov-Smirnov test to evaluate the significance of 
the differences between the ^-distributions of binary and (apparently) non-binary objects, 
finding a 4.7% probability that the two populations were drawn from the same distribution. 
Here we redo the Kolmogorov-Smir nov test, taking i nto account the time-variability of the 
inclinations. Of the 101 objects in iNoll et al.l ( l2008l ). 93 are in our sample of CKBOs; 21 
of these are binary. We run the K-S test 1000 times, each time randomly assigning values 
of the inclinations from a normal distribution having average and standard deviation in i 
found in our 10 Myr simulation. This calculation finds that there is a 6.5 ± 3% probability 
that the binary and non-binary CKBOs are derived frorn the same overall inclination dis- 
tribution. This result is consistent with iNoU et al.l (120081 ) 's results. We conclude that the 
time-variability of inclinations in the classical Kuiper belt does not have a strong effect on 
the inclination dependency of the binary fraction. 



5. Summary 

We studied the dynamics of the presently observed classical and resonant Kuiper belt 
objects, with a focus on their inclination evolution. 

1. Planetary perturbations over gigayear timescales can produce inclination changes in 
the orbits of Kuiper belt objects. We find that large changes {Ai > 5°) are confined 
to a small fraction of the parameter space, and they typically last only a relatively 
small fraction of the time over ~ 4 gigayears. Most objects (90% of CKBOs and 77% 
of resonant KBOs) remain within 5° of their current inclinations and have short term 
variations with amplitude typically < 2° (relative to the invariable plane). Overall, 
the debiased inclination distribution of the non-resonant classical Kuiper belt is nearly 
stationary under planetary perturbations over gigayear timescales. 

2. Considering the inclination distribution of classical Kuiper belt objects, we find that 
there is some transfer between the high- and low-z populations, but at any given time 
relatively few CKBOs (< 5%) will be found at inclinations inconsistent with their 
original inclination classification. Most CKBOs that experience changes in i larger than 
5° are ejected from the classical belt by close encounters with Neptune, but some objects 
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(mostly in the outer half of the classical belt) can undergo large inclination changes 
while maintaining fairly low eccentricities thus avoiding close Neptune encounters. 

3. We identify many KBOs within the classical belt region that are in high-order MMRs 
with Neptune. We find that some current members of the 7:4 and 9:5 MMRs escape 
from these resonances into the non-resonant classical belt population without having 
destabilizing close encounters with Neptune. Many of the CKBOs in our simulation 
that experience large changes in i are located near (but not in) these high-order reso- 
nances. Temporary capture into high-order MMRs could explain the CKBOs with the 
most changeable inclinations. 

4. The resonant KBOs are more likely to undergo large changes [Ai > 5°) in inclination 
than non-resonant objects. The largest variations in i can be attributed to the Kozai 
mechanism operating in conjunction with a mean motion resonance. We estimate 
that, at any given epoch, 5 to 10% of observed resonant objects will be at inclinations 
inconsistent with their original inclination classification. 

5. We reassessed the previously proposed correlations between inclination and the bina- 
rity, sizes, and colors of classical Kuiper belt objects, using updated observations and 
taking account of the time- variability of the inclinations. 



The inclination time- variability does not greatly affect iNoll et al.l ( l2008l )'s finding 
that the binary and non-binary CKBOs have different inclination distributions; we 
find only a 6.5±3% probability that the binary and non-binary CKBOs are derived 
fro m the same inclination distribution, very similar to the ~ 5% probability found 



bv lNoU et all tom . 

(b) For the observationally near-complete population of bright objects (absolute mag- 
nitudes H < 5), we find a 4 — 7% probability that the large objects have the same 
inclination distribution as the rest of the CKBOs. This is a slig htly higher prob 



abilit y than the previously reported 3% chance of similarity (ILevison fc Stern 



20011 ). although the two results are consistent within our error estimates for the 
probability; most of the difference is due to the larger data set available since the 
previous study. Using the Spearman rank test, we find no significant correlation 
between i and H for the CKBOs as a whole or for the bright {H < 5) CKBOs 
considered separately. 

We find a decreased significance of the correlation between color (represented 
here by spectral gradient) and inclination of CKBOs; previous estimates reported 
the correlation coefficient Rg = —(0.6 — 0.7) (8cr significance) while our analysis 
finds Rg = —0.46 ± 0.07 (4 ± 0.5a significance). A large part of the difference 
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between our result and the previous results is due to the larger data set of observed 
objects, with a small (but significant) contribution from taking account of the 
observational uncertainty in the measured colors and the time-variability of i. 
There is also evidence that the Haumea coUisional family is responsible for at least 
part of the apparent correlation; removing this high-inclination, neutrally colored 
group of objects lowers the correlation and its significance to Rg = —0.36 ± 0.07 
(3.1 ± 0.6a significance). 

This research was supported by grant no. NNX08AQ65G from NASA's Outer Planets 
Research program. 
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Table 1. Resonant Objects 



MMR 



Designations of Members 



3:2 



7:4 



2:1 



5:3 



9:5 



1993 RO 


1993 SB 


1993 SC 


1994 JRi 


1994 TB 


1995 HM5 


1995 OYqfr) 


1995 QZ9 


1996 RR20 


1996 SZ4 


1996 TPrfi 


1996 TOrfi 


1997 OT/i 


1998 HHi^i 


1998 HKi.;i 

tj tj \J XXXVJ^JdJ^ 


1998 HOiri 


1998 URac! 


1998 USac! fr") 


1998 YGaa 


1998 WScii 


1998 WTJqi 


1998 WVqi 


1998 WWo/i 

-Ltyt/u vy yv ^4 


1998 WZqi 


1999 CEiio 


1999 CMi fr") 


1999 RKois 


1999 TCqk 


1999 TR11 


2000 CKinr 


2000 EBi7Q 


2000 FBs 


2000 FVsQ {t\ 


2000 GEi/17 


2000 GNi7i 


2000 YHn 


2001 FTjin/i 


2001 FRiQK 


2001 FTI179 


2001 KB77 


2001 KDt7 


2001 KN77 fr") 


2001 KO77 


2001 KX7fi 


2001 KY7K 


2001 0F9QS 


2001 OG9QS 


2001 0H9no 


2001 RUiA-i 


2001 RXiA'? 


2001 UOie frl 


2001 VN71 


2001 YTi^n 


2002 GE9r;i fr") 


2002 CW99^ 


2002 GE w 


2002 GF19 


2002 GLc!9 


2002 GV-!9 


2002 GW'^i 


2002 GY32 


2002 VDicis 


2002 VE95 


2002 VRi9s 


2002 VU130 


2002 VX130 


2002 XV93 


2003 AZ84 


2003 FBi9s 


2003 FFi9s 


2003 FLn7 


2003 HAt;7 


2003 HDt;7 


2003 HFr;7 


2003 OBni 


2003 QH91 


2003 QXiii 


2003 SO317 


2003 SR317 


2003 TH.;s fr) 


2003 UTooT 


2003 UVoQo 


2003 UZ/iiq 


2003 VS9 


2003 WAioi 


2003 WTJi^T fr") 


2004 DW 


2004 EHnc 


2004 ElnK 


2004 FVnr 


2004 EWq^ iv) 


2004 FUiAs 


2004 FW1RZ1 


2004 USin 


2004 VT75 


2004 VZtk 


2005 EZoQfi 


2005 EZqnn ir) 


2005 GAi«7 


2005 GBi«7 


2005 GEis7 


2005 GFis7 


2005 GV210 


2005 TV18Q 


2006 HJ123 


2007 .TFn 


2007 .THic fr^ 


Pluto 






(1994 GVg) 


1999 CDi.;s fr) 


1999 CO153 


1999 HG12 


(1999 HRii) 


1999 HTii 


1999 KR18 


1999 RH215 


2000 FD8 (r) 


(2000 FX53) 


2000 OP67 


2000 OY51 (r) 


2000 YUi (r) 


2001 HA59 (r) 


2001 KJ76 


2001 KO76 


2001 KP76 (r) 


2001 KP77 (r) 


2001 QE298 


2002 PA149 (r) 


2002 PB171 


2003 QWiii 


2003 QX91 


(2003 YJ179) 


2004 OK14 (r) 


2004 PW107 (r) 


2004 VU75 (r) 


2005 SF278 (r) 






1996TR66 (r) 


1997SZio (r) 


1998SMi65 (r) 


I999RB215 (r) 


I999RB216 (r) 


2000JG81 (r) 


2OOOQL251 


2OOIFQ185 


2OOIUP18 


2002PUi7o (r) 


2OO2WC19 


2002VDi3o 


2003FEi28 (r) 


2OO4TV357 


2OO5CA79 


2OO5RS43 










1994 JS (r) 


1997 CV29 


1999 CX131 (r) 


2000 PL30 (r) 


2000 QN251 (r) 


2001 XP254 (r) 


2001 YHi4o (r) 


2002 GS32 


2002 VA131 


2002 VV130 


2003 US292 (r) 


2003 YW179 


2005 SE278 


2005 TN74 (r) 


2006 QQ180 


1999 CS153 (r) 


2000 CN105 (r) 


(2000 FR53) 


2000 YZi (r) 


2001 KL76 


2002 GD32 


(2004 VS75) 


(2005 JH177) 
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Table 1 — Continued 



MMR Designations of Members 



4:3 


1995 DA2 
2004 TX357 (r) 


1998 UU43 
2005 ER318 (r) 


1999 RW215 


2002 FWe 2003 SS317 (r) 


5:4 


1999 CP133 
2005 SC278 


2001 XH255 


2002 GW32 (r) 


2003 FC128 2003 QB92 


11:6 


1999 CP153 (r) 


2001 KU76 


2005 EB318 (r) 




8:5 


2005 VZ122 (r) 








11:8 


2001 XS254 (r) 








16:9 


2005 GCi87 (r) 









Not e. — Entries in bold f ace in dicat e objects not previou sly identified as resonant bv lGladman et al 



( 2008h . Lvkawka fc Mukai ( 2007 ). or Chiang et al. ( 2003 ). Entries in parentheses are objects whose 
test particles exit the resonance by the end of the 4Gyr simulation without being removed from the 
classical belt region. Entries followed by (r) indicate objects whose test particles are removed from 
the simulation after close encounters with Neptune. 
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Fig. 1. — Locations of identified mean motion resonances in the classical Kuiper belt. The 
vertical height associated with each resonance is arbitrarily scaled to the order of the res- 
onance. The numbers above the lines are how many objects we have identified in each 
resonance. 
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Table 2. Kuiper Belt Surveys 



Number 


Survey 


Reference (s) 


of Objects 






171 


Deep Ecliptic Survey 


Elliot et al. (2005) 


34 


Canada-France Ecliptic Plane Survey 


Jones et al. (2006) 






Kavelaars et al. (20091 


33 


Canada-France-Hawaii Telescope Survey 


Truiillo et al. (2001a) 


8 


Spacewatch 


Larsen et al. (2001) 






Larsen et al. (2007) 


7 


Mauna Kea 


Jewitt et al. ( 1998) 


6 


Distant Disk Survey 


Allen et al. (2002) 


19 


miscellaneous 


Irwin et al. ( 1995) 



Jewitt fc Luu (1995) 
Jewitt et al. (1996) 
Gladman et al. (1998) 
Luu fc Jewitt (1998) 
Truiillo fc Jewitt (1998) 
Gladman et al. (2001) 
Truiillo et al. (2001b) 
Truiillo fc Brown (2003) 



Note. — Surveys used in this study to construct a debiased model of the non- 
resonant classical Kuiper belt. 



- 30 - 



0.30 



0.25 



0.20 



0.15 



0.10 



0.05 



DES CKBOs 

this work (using the DES definition of CKBOs) 



10 15 
inclination (deg) 



20 



25 



Fig. 2. — Comparison of our debiasing procedure to that from the Deep Echptic Survey 
(DES) JCulbis et ahlboiol ). The distributions contain 150 DES KBOs that meet the DES def- 



Gulbis et al. 



Gulbis et al. 



inition of CKBO and are present in the data sets for this work as well as that of 
( 2010 ). The dark line shows the debiased ^-distribution for these KBOs from 
(120101 ) and the grey line shows our debiased distribution for the same objects. The debiased 
z-distribution resulting from the definition of CKBOs used in this work is shown in Figure [71 
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Fig. 3. — Variation in inclination vs. average inclination for the initial 10 Myr integration 
relative to the ecliptic plane (black) and the invariable reference plane (gray). 
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Fig. 4. — Variation in inclination as a function of average semimajor axis for the non-resonant 
(black) and resonant (gray) test particles over the 10 Myr integration. The left panel shows 
test particles that survive the 4 Gyr simulation without a close encounter with Neptune, and 
the right panel shows those that are removed from the simulation by a close encounter. 
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Fig. 5. — Variation in inclination as a function of average semimajor axis for the non-resonant 
(black) and resonant (gray) test particles over the 4Gyr integration. The left panel shows 
test particles that survive the 4 Gyr simulation without a close encounter with Neptune, and 
the right panel shows those that are removed from the simulation by a close encounter. The 
50 Myr time period prior to particle removal is excluded from the calculation of the average 
and rms variation. 
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Fig. 6. — Maximum range in inclination as a function of average semimajor axis for the 
non-resonant (black) and resonant (gray) test particles over the 4 Gyr integration. The left 
panel shows test particles that survive the 4 Gyr simulation without a close encounter with 
Neptune, and the right panel shows those that are removed from the simulation after a close 
encounter. The 50 Myr time period prior to particle removal is excluded from the calculation 
of the average and maximum range. 
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Fig. 7. — Snapshots of the debiased inchnation distribution for the non-resonant CKBOs 
taken at the beginning of the simulation, at 500 Myr, and at 4 Gyr. (Note that the i- 
distribution depends heavily on the definition used for CKBOs. The differences between this 
figure and Figure p are entirely owed to the differences between our classification scheme 
and that used by (iGulbis et al.l 120101 ).) The two smooth curves represent model fits of the 
CKBO inclination distribution to the sum of two Gaussians multiplied by sini (Equation |3]). 
The solid grey line is our bes t-fit model to the initial distribution, and the dashed grey line 
is the model of iBrownl (120011 ) for comparison. 
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Fig. 8. — Time- averaged fraction of the total observed resonant population (dashed line) and 
of the total debiased non-resonant population (solid line) that is misclassified. This fraction 
is plotted as a function of the cutoff between the low- and high-inclination populations. Left 
panel: inclinations referred to the invariable plane; Right panel: inclinations referred to the 
ecliptic plane. 
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Fig. 9. — Eccentricity, inclination, semimajor axis, and perihelion distance evolution over 4 
Gyr for a test particle that maintains a low i throughout the entire simulation. 
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Fig. 10. — Eccentricity, inclination, semimajor axis, and perihelion distance evolution over 
the first 750 Myr of the 4 Gyr simulation for a test particle that has a distant encounter 
with Neptune (at about 500 Myr) and a consequent jump in a, e, and i. 




Fig. 11. — Eccentricity, inclination, semimajor axis, and perihelion distance evolution over 
4 Gyr for a test particle with an anomalously long-period variation in i and e. 
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Fig. 12. — Eccentricity, inclination, resonance angle {(f) = 7 Xkbo — — SiUkbo) , and argument 
of perihelion, u, evolution over 4 Gyr for a test particle exhibiting intermittent libration 
within the 7:4 MMR with Neptune, as well as intermittent oj libration. 
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Fig. 13. — Spectral gradient (percent reddening per 100 nm) vs. average inclination over the 
10 Myr simulation. The error bars for the inclination are the rms variation in i over 10 Myr. 
The 6 objects enclosed in the ellipse are all members of the Haumea collisional family. 



